{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generateGalaxy():#生成两个星系，每个星系包含100个恒星，假设恒星只与两恒星中心黑洞有万有引力作用，不考虑恒星间作用\n",
    "    w0=np.zeros(808) #定义一个0数组，用来存储各个恒星初始值，0-399为星系A，400-799为星系B,800-803记录星系A中心,803-807记录星系B中心\n",
    "    \n",
    "    w0[800]=-2\n",
    "    w0[801]=0.15\n",
    "    w0[802]=3  #星系A的y轴\n",
    "    w0[804]=2\n",
    "    w0[805]=-0.15\n",
    "    w0[806]=-3  #星系B的y轴\n",
    "    r=np.linspace(0.5,2,4)\n",
    "    num=0\n",
    "    for i in r: #星系A\n",
    "        thetaR=np.arange(0,i*2*np.pi,0.1*np.pi)\n",
    "        theta=thetaR/i\n",
    "        v=np.sqrt(0.1/i) #能保持圆周运动的速度，0.1是系数k，若该则在对于的传入odeint函数里改\n",
    "        for j in theta:\n",
    "            w0[num]=i*np.sin(j)+w0[800]    #分别是一个恒星的x,vx,y,vy\n",
    "            w0[num+1]=v*np.sin(np.pi/2+j)+w0[801] #这里+0.1是为了个星系A，星系B一个初始选转\n",
    "            w0[num+2]=i*np.cos(j)+ w0[802]\n",
    "            w0[num+3]=v*np.cos(np.pi/2+j)\n",
    "            num=num+4\n",
    "\n",
    "    for i in r:  #星系B\n",
    "        thetaR=np.arange(0,i*2*np.pi,0.1*np.pi)\n",
    "        theta=thetaR/i\n",
    "        v=np.sqrt(0.1/i) \n",
    "        for j in theta:\n",
    "            w0[num]=i*np.sin(j)+w0[804]    #分别是一个恒星的x,vx,y,vy\n",
    "            w0[num+1]=v*np.sin(np.pi/2+j)+w0[805]\n",
    "            w0[num+2]=i*np.cos(j)+w0[806]   #星系A中心的初始坐标为（0，0），这里把星系B的坐标设为（0，7）\n",
    "            w0[num+3]=v*np.cos(np.pi/2+j)\n",
    "            num=num+4\n",
    "\n",
    "    return w0\n",
    "    \n",
    "def galaxyOrbit1(w,t):\n",
    "    temp=[]\n",
    "    ax=w[800]\n",
    "    vax=w[801]\n",
    "    ay=w[802]\n",
    "    vay=w[803]\n",
    "    \n",
    "    bx=w[804]\n",
    "    vbx=w[805]\n",
    "    by=w[806]\n",
    "    vby=w[807]\n",
    "    \n",
    "    rab=np.sqrt((ax-bx)**2+(ay-by)**2)\n",
    "    cosAB=(ax-bx)/rab\n",
    "    sinAB=(ay-by)/rab\n",
    "    \n",
    "    \n",
    "    for i in range(0,200):\n",
    "        j=4*i\n",
    "        k=0.1 #计算加速度系数\n",
    "        \n",
    "        ra=np.sqrt((w[j]-ax)**2+(w[j+2]-ay)**2)\n",
    "        cosA=(ax-w[j])/ra\n",
    "        sinA=(ay-w[j+2])/ra\n",
    "        rb=np.sqrt((bx-w[j])**2+(by-w[j+2])**2)\n",
    "        cosB=(bx-w[j])/rb\n",
    "        sinB=(by-w[j+2])/rb\n",
    "\n",
    "        dx=w[j+1]\n",
    "        dvx=k*cosA/ra**2+k*cosB/rb**2 #恒星被A,B中心黑洞吸引\n",
    "        dy=w[j+3]\n",
    "        dvy=k*sinA/ra**2+k*sinB/rb**2     \n",
    "        \n",
    "\n",
    "        temp.append(dx)\n",
    "        temp.append(dvx)\n",
    "        temp.append(dy)\n",
    "        temp.append(dvy) \n",
    "    \n",
    "    dax=vax\n",
    "    dvax=-k*cosAB/rab**2\n",
    "    day=vay\n",
    "    dvay=-k*sinAB/rab**2\n",
    "    \n",
    "    dbx=vbx\n",
    "    dvbx=k*cosAB/rab**2\n",
    "    dby=vby\n",
    "    dvby=k*sinAB/rab**2  \n",
    "    \n",
    "    temp.append(dax)\n",
    "    temp.append(dvax)\n",
    "    temp.append(day)\n",
    "    temp.append(dvay)\n",
    "\n",
    "    temp.append(dbx)\n",
    "    temp.append(dvbx)\n",
    "    temp.append(dby)\n",
    "    temp.append(dvby)\n",
    "    return temp\n",
    "\n",
    "def plotGalaxy(w):\n",
    "    rows,col=w.shape\n",
    "    for row in range(rows):\n",
    "        plt.clf()\n",
    "        for i in range(100):\n",
    "            j=i*4\n",
    "            plt.plot(w[row][j], w[row][j+2], \"ro\")\n",
    "        \n",
    "        for i in range(100,200):\n",
    "            j=i*4\n",
    "            plt.plot(w[row][j], w[row][j+2], \"bo\")\n",
    "        plt.plot(w[row][800], w[row][800+2], \"go\")\n",
    "        plt.plot(w[row][804], w[row][804+2], \"go\")\n",
    "        plt.axis(\"equal\") #xy轴间距相等\n",
    "        plt.xlim((-16,16))\n",
    "        plt.ylim((-16,16))\n",
    "        plt.savefig(\"{:0>3d}\".format(row)+\".png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHllJREFUeJzt3XuQXGeZ3/Hvb0YWMLLBWBIGjNXjLM5WGRy8eGJCOUuZNRijgvWSgi0rg6MFKrOSocrUJqmFnRQmbKZqN8mGGBabGsCLgppbNuvFBHGxnSWEKm4jYrDMJQh7Rgi78FgGG6ENtqQnf5xuT3fPOd2nL6dv8/tUdXX36bdPv60pnaff2/MqIjAzM6uaGHQFzMxsuDgwmJlZHQcGMzOr48BgZmZ1HBjMzKyOA4OZmdVxYDAzszoODGZmVseBwczM6mwadAU6sW3btpienh50NczMRsrBgwcfjojtrcr1JDBIuhV4DfBQRLywcuwc4FPANLAM/H5E/DzlvbuBf1t5+u8jYl+rz5uenmZpaakXVTcz2zAkreQp16uupI8CVzccewdwV0RcCNxVeV6nEjxuBF4CXAbcKOmZPaqTmZl1oCeBISK+AjzScPgaoPrrfx/weylvfRVwR0Q8UmlN3MH6AGNmZn1U5ODzuRHxIEDl/lkpZc4DflLz/Gjl2DqS5iQtSVpaXV3teWXNzCwx6FlJSjmWmgc8IhYjYiYiZrZvbzl2YmZmHSoyMPxM0nMAKvcPpZQ5Cpxf8/x5wAMF1snMzFooMjDcDuyuPN4NfCalzBeBqyQ9szLofFXlmJmZDUhPAoOkTwBfA35T0lFJbwH+DHilpB8Br6w8R9KMpA8DRMQjwJ8C36rc3lM5ZmZmA6JR3NpzZmYmvI7BzKw9kg5GxEyrcoMefDYzsyHjwGBmZnUcGPqgXIbpaZiYSO7L5UHXyMwsmwNDhmYX83Yu9OUyzM3BygpEJPdzcw4OZja8PPiconoxP3Fi7djUFCwuJo+zXpudXX+u6ekkGDQqlWB5uZe1NjNrLu/gswNDimYXc2jvQj8xkbQUGklw+nQ3tTQza49nJXXhyJHs481eS7NjR3vHzcwGbcMGhmbjBM0u5u1e6BcWkq6mWlNTyXEzs2G0IQNDqwHhZhfzdi/0s7PJ+EOplHQflUrZ4xFmZsNgQ44x5BkQLpdhfj7pItqxI7nwVy/mzV4zMxtWHnxuoqgBYQcMMxtmeQNDT/Z8HiXlchIYTp1a/1o3A8KNU1yr3VPg4GBmo2VDjTFUL95pQaHbAeH5+fq1DZA8n5/v/JxmZoOwoQJD2sUbYHKy+wHhdqexmpkNqw0VGLIu0qdP5wsKnUxxPeec9vMkObeSmQ3S2AeG2ovsRMa3zTO20MkUV4Bjx9rLk+TcSmY2aGMdGBovst2MLbQaQ6hdrwDJDKc0zcYdymXYvbuNsQo3LcysCBFR2A34TeDumttjwNsbylwBPFpT5l2tznvppZdGHqVSRBIS6m+TkxFS8vr+/blOFVL6uaT8n1t7a/zc/fsjpqayy6/7nLQ3TE3l/0JmtuEAS5Hj2t23dQySJoGfAi+JiJWa41cA/zoiXpP3XHnXMfRyvUKzRXELC8kv+pWVZCA7rWXS6Iwz4OlPT7qa8rxnXZI+p201szYNYxK9K4Ef1waFom3Zkn78nHPaP1faGIKUXJuvu27tGp0nKAA88UQSFPK8p667q9p9lBYUwNOgzKxr/QwM1wKfyHjtpZK+I+nzkl7Qiw8rl+H48V6cKZE2hlBtjRTZ6KqbSls7aJLFaVvNrEt9CQySNgO/C/y3lJe/DZQi4kXA+4G/zTjHnKQlSUurq6stP7PZwrJHHslR6Rb60QMnwb59yePpaZh44y6mT9zL9byfae5nglNMcz9ldiWFnLbVzHqgL2MMkq4B3hoRV+UouwzMRMTDWWXyjDFkjS9AZ93wabu69Utt6yQRwNq0pyl+xeLWdzJ700ucf8PMMg3bGMMuMrqRJD1bSiZ3SrqsUqdj3X5gVo+K1N6P6mqX/hvfOJigAGkBrn4u7Am2MH/m+3IHBc9yNbNmCg8MkqaAVwJ/U3Nsj6Q9laevBw5J+g7wPuDa6EEzJmuweM+e+utns4tkuQxvfnPzLv1hkXfM2QvozKyVsU673SoNdlr30NTU2mDvtm1rM4eGXd7uMc9yNdu4vB9DDq0uklmrl4dNbTBrpai9KMxs+A3bGMNQ6iYj6rAEDSlJo5F3zLndPavNbOPZ0IGh2UXy+uuz37d+ltDgRMCBA/nLt7tntZltPBs6MGRdJJ//fLjllvT3nHHG8ASFqnYWO9cu1JOS+273ojCz8bKhA0PWRfLLX85+z1/9FWzdmvMDLi7D26fhxonk/uJ8U3+q6cFLpXyf1W430OxsMoZy+nRy76BgZrU2dGCA9Itkq9xFjz2W48QXl+G1c3D2CiiS+9fOpQaHxvGK06fXunduuil9n4eqzTzOws6v5qiQmVk+Gz4wpJmczD4+P58kwGvpynnY3LAibvOJ5HiDtK6p6h4M1VZNVsvhFIIPf8QLEcysZxwYUlxxRfrxubk2+vOfkVEw63iK6mfNzsKZZ6aXOcUZzD9xY/PkUGZmbXBgaFAuw9e+tv74lVfCzTe30Z//aEbBlON5thw9spI94n2EHU63bWY948DQIG0LT4DDh5P7rL2d17lrAR5vKPj4VHK8QdrCstoppOXrv8oE2QMfOzjihQhm1jMODA1aLXqbnU0WlLVc4HbPLHx2EX5RglBy/9nF5HgLW7euTSEtl2Hugy/mFJtSy27m1yyc8e+8EMHMemZDp8RIkyeXULMN1Hoh32cFW3mYm7b+qdNtm1kuTonRoayuouPH1yb+FN2dX3v+rM8SwcOlf8zsw/nTbZuZ5eHAUFFNv33ddfC0p63fL/rYsbX01EV35+/YsVafrAbdDh1195GZFcKBgfV7FBw7lj4AXV1bkHsAugNTU7BzZ/Otnad0goU9R9xSMLNCODCQPhMp65f6kSP1qTSy1E5BbScT6+JikhQva7e4UgkWPzbF7M3/NP9Jzcza4MBAe2MGEUkXDzTvyamdgpp3fH9yMunKymopSM5tZGbFS58D2UOSloFfAqeAk40j4pX9nm8CdgIngD+IiG8XXa9aO3akX4yz0muvrMCb3tT7LKutcjR5qYKZ9UO/Wgwvj4hLMqZJvRq4sHKbAzISXhcnK/32nj3Z3UVPPAEnTxZft9r6eKzZzPphGLqSrgH+ayS+Dpwt6Tn9rEBW+u3LL+9nLdJ5zwQz67d+BIYAviTpoKS5lNfPA35S8/xo5VhfNabfhuYzg/qhVPKeCWbWf/0IDJdHxItJuozeKullDa+nzdlZ13svaU7SkqSl1dXVIupZJytnUr+kdh1VFzdMTCT3TrVtZgUoPDBExAOV+4eA24DLGoocBc6vef484IGU8yxGxExEzGzfvr2o6j6p2UylrP0aurF5c5IjKbPrqHGxxcrK2oo7M7MeKjQwSNoi6azqY+Aq4FBDsduBf6HEPwEejYgHi6xXHlkzgEol2Lev+QK3pz619ZacW7bUB4Jbb4WHH27SdZTWhKmuuDMz66Gip6ueC9yWzEhlE/DxiPiCpD0AEfFB4ADJVNXDJNNV31RwnXJZWEh+kNdei6vdO9WL9vx8/RjE5GTynptvTp5PTGRPaT1+vM0KtUr7ambWI86u2kS5nFz8j1S2O6gNCnnkydQ6mJOZ2Ubk7Ko90DhTqd2ZQVnrIzpaj9DTk5mZZXNg6LWamUOz89Ms7v7quvURHU09zVps4XmsZtZj7krqperMocaBCV/AzWwIuCtpEDxzyMzGgANDLw3RzCGvhTOzTjkw9FLW4oc+p0X1Wjgz64YDQy8Nycwh92iZWTccGHppSGYODVGPlpmNoMI36tlwZmcHPgMpa+Mhb/RjZnm4xTCGhqRHy8xGlANDn/VjttCQ9GiZ2YhyV1IfNa5/q84Wgt5ftIegR8vMRpRbDH3k2UJmNgocGPrIs4XMbBQ4MPTRkKx/MzNryoGhjzxbyMxGgQNDH3m2kJmNgsICg6TzJf2dpO9LulfSDSllrpD0qKS7K7d3FVWfYdHt5j9mZkUrcrrqSeBfRcS3JZ0FHJR0R0R8r6Hc/46I1xRYDzMza0NhLYaIeDAivl15/Evg+8B5RX2emZn1Rl/GGCRNA78FfCPl5ZdK+o6kz0t6QT/qY2Zm2Qpf+SzpTOC/A2+PiMcaXv42UIqI45J2An8LXJhxnjlgDmCH53eamRWm0BaDpDNIgkI5Iv6m8fWIeCwijlceHwDOkLQt7VwRsRgRMxExs3379iKrbWa2oRU5K0nAR4DvR8R/zijz7Eo5JF1Wqc+xouo0Lrxtp5kVqciupMuB64B7JN1dOfYnwA6AiPgg8Hpgr6STwN8D10ZEFFinkdfPRHxmtjFpFK/DMzMzsbS0NOhqDMT0dPomPKVSsi7CzCyLpIMRMdOqnFc+jxgn4jOzojkwjBgn4jOzojkwjBgn4jOzojkwjBgn4jOzonlrzxHkbTvNrEhuMVhveZGF2chzi8F6x4sszMaCWwzWO/Pza0Gh6sSJ5LiZjQwHBusdL7IwGwsODNY7XmRhNhYcGKx3vMjCbCw4MFjveJGF2VjwrCTrLS+yMBt5bjGYmVkdBwYzM6vjwGBmZnUcGMzMrE7hgUHS1ZJ+KOmwpHekvP4USZ+qvP4NSdNF18nMzLIVGhgkTQIfAF4NXATsknRRQ7G3AD+PiOcD7wX+vMg6mZlZc0W3GC4DDkfEfRHxOPBJ4JqGMtcA+yqP/xq4UpIKrpeZmWUoOjCcB/yk5vnRyrHUMhFxEngU2FpwvczMLEPRgSHtl390UAZJc5KWJC2trq72pHJmZrZe0YHhKHB+zfPnAQ9klZG0CXgG8EjjiSJiMSJmImJm+/btBVXXzMyKDgzfAi6UdIGkzcC1wO0NZW4Hdlcevx74nxGxrsVgZmb9UWiupIg4KeltwBeBSeDWiLhX0nuApYi4HfgI8DFJh0laCtcWWSczM2uu8CR6EXEAONBw7F01j/8f8Iai62FmZvl45bOZmdVxYDAzszoODGZmVseBwczM6jgwmJlZHQcGM+tYuQzT0zAxkdyXy4OukfWC93w2s46UyzA3BydOJM9XVpLn4G2/R51bDGbWkfn5taBQdeJEctxGmwODmXXkyJH2jtvocGAws47s2NHecRsdDgxm1pGFBZiaqj82NZUct9HmwGBmHZmdhcVFKJVASu4XFz3wPA48K8nMOjY760AwjtxiMLPx5EUWHXOLwczGjxdZdMUtBjMbP15k0RUHBjMbP15k0RUHBjMbP15k0ZVCAoOk/yjpB5K+K+k2SWdnlFuWdI+kuyUtFVEXM9uAvMiiK0W1GO4AXhgR/wj4v8A7m5R9eURcEhEzBdXFzDYaL7LoSiGzkiLiSzVPvw68vojPMTPL5EUWHevHGMObgc9nvBbAlyQdlDTX7CSS5iQtSVpaXV3teSXNzCzRcYtB0p3As1Nemo+Iz1TKzAMngayVJZdHxAOSngXcIekHEfGVtIIRsQgsAszMzESn9TYzs+Y6DgwR8Ypmr0vaDbwGuDIiUi/kEfFA5f4hSbcBlwGpgcHMzPqjqFlJVwN/DPxuRJzIKLNF0lnVx8BVwKEi6mNmZvkVNcbwl8BZJN1Dd0v6IICk50o6UClzLvBVSd8Bvgl8LiK+UFB9zMwsp0ICQ0Q8PyLOr0xDvSQi9lSOPxAROyuP74uIF1VuL4gITzA2s+GyQRPxOYmemVmaDZyIzykxzMzSbOBEfA4MZmZpNnAiPgcGM7M0GzgRnwODmVmaDZyIz4HBzCzNBk7E58BgZqOpH1NJZ2dheRlOn07uN0BQAE9XNbNRtIGnkvaDWwxmNnqGaCrpOK6Bc4vBzEbPkEwlHdeGi1sMZjZ6hmQq6RA1XHrKgcHMRs+QTCUdkoZLzzkwmNnoGZKppEPScOk5BwYzG01DMJV0SBouPefAYGaWV8MUpFnKw9Bw6TkHBjMbaz2bTlqdgrSyAhFPTkGapTzohkvPOTCY2djKuJZ3Fhx6NQVpBBY+FBYYJL1b0k8rW3veLWlnRrmrJf1Q0mFJ7yiqPma28fR0OunKSvrx2ilI118PmzYl/UqbNiXPa/U0UhWn6BbDe2u29zzQ+KKkSeADwKuBi4Bdki4quE5mtkF0PJ209lf9tm1w5pnZZSOSMk95CtxyC5w6lRw/dSp5ftZZaxf+EVn4MOiVz5cBhyPiPgBJnwSuAb430FqZ2VjYsSP9h37T6aSNy5mPHWv9Qc3KHD++thx6RBY+FN1ieJuk70q6VdIzU14/D/hJzfOjlWPrSJqTtCRpaXV1tYi6mtmYyT2dtLaFsHv3+l/13TpxIjlvRPrrQ7bwoavAIOlOSYdSbtcAtwC/AVwCPAj8RdopUo6l/stFxGJEzETEzPbt27uptpltELnWwTX2+1e7gnot67xDuPChq8AQEa+IiBem3D4TET+LiFMRcRr4EEm3UaOjwPk1z58HPNBNnczMaqWugyu6hZBXbaQaotlKhY0xSHpORDxYefo64FBKsW8BF0q6APgpcC3wz4uqk5nZujGEoloI7RiyNK2KrD6vbk8sfYykGymAZeAPI+JBSc8FPhwROyvldgL/BZgEbo2Ilm2qmZmZWFpaKqTeZjbmpqezp54OwtQUPO1p6QPYpVLSzOkRSQcjYqZluaICQ5EcGMysYxMT2YPARdi7F/bt66y7Skr6wHokb2Dwymcz2xiqffhZQWFyspjPPXAgGccoldp/74BmKzkwmNn4K5cpv+lOple+zASnmOZ+yuxae31qKunTb5zb2gsrK0mLYWEhOzhs3TpUaVodGMxs7JVv+AZzT/wlK0wTTLDCNHN8KAkO1ZlBBw4UNzupuro5bWFFVbVVMQRpWj3GYGZjb1rLrDC97niJZZajcrzosYfqeEG5DDfcsH6weWqq8GDgMQYz29DKZZjedpwJnWaF9C6cI9T04Rfdn189/+xseu6lIcqZ5MBgZmOnXIa5N59k5diZBBOkJ1mAHVtruo527kx+1bc698Uw/XaYuDG5L1+co0JScv6qIc+Z5MBgZmNnfh5OPN58/e7Upl+zcFPll3u5nAwQt+hKKl8Mc6+FlbMhlNzPvTZHcIhIzl9dzTzkm0U7MJjZ2Mn+4R2I05RYZvEZ/2atOz8tHXaK+SvhxOb6Yyc2J8dbqu0qSmudDFHOpEGn3TYz67msdNslVljmguTJIwLelzzO2YVz5BntHV9f8EjSavjwh9e3TnbvHpp9Qd1iMLOxs7AAU5tP1h2b4lcs8CdrB845Z+1xzi6cHY+2d3x9wR3JjKQnnlj/2qc/nfMkxXNgMLOxMzsLi2/5OiVW1rqO+JfM8om1Qo89ttbn32x9QY2Fu2Dq8fpjU48nx1s644zkc7I29cmzIVCfODCY2ViaPfBGlpnmNJMsc0F9UIDkV/sNNyRpMq67LlnH0Oqc98DiZ6H0C1Ak94ufTY639PSnt+4qGpLU217gZmbDqVxOBmuPHEm6YBYW2uuD73eyvFYk2LMn2Qc6zZYtSX1rB8F7vOjNC9zMbHQ17qpW3Z+gnV/QQzL180lTU9lBYfNmeOpT18+MGtCiNwcGMxs+adNH271I5lywlku3mVcl+NWvsl+/9VZ45JH01waw6M2BwcyGT7crg3MuWMut213eWtVjdnaoFr05MJjZ8On2IllpcZTZxTT3p6faHhbV1kjazKi0RW/9GKCOiJ7fgE8Bd1duy8DdGeWWgXsq5Zbynv/SSy8NMxtj+/dHTE1FJL+1k9vUVHI8Dyn2syumOF5/Co7HfnbVn3fQt7176793qRQhJfeN33fv3uS1Dv9d8l5nCwkMdR8AfwG8K+O1ZWBbu+d0YDDbAFpdJJsplaLE/anX4RL3Dz4YQMTkZH1QyPPv0RgUnvxSpVynyBsYCk2JIUnA7wO/U+TnmNkYmp3taJpmuQzzxw+xwpbU159MtS0NZjprp1NQ5+ez69vjAeqixxh+G/hZRPwo4/UAviTpoKS5gutiZmPuyVmux84kM9U2R5Id0j72Mdi/v3czl/LoZme2Zhf/Hg9QdxwYJN0p6VDK7ZqaYrugcblhncsj4sXAq4G3SnpZk8+bk7QkaWl1dbXTapvZGGuVJHVqChb2T8PycnJg9+7+tBp6EXxqczs16nFW1sJWPkvaBPwUuDQijuYo/27geET8p1ZlvfLZzNI0W+xcKtUsnq42LYra4xnWuqoau6w67Urati09n9KZZ8Ivf5mzSoNf+fwK4AdZQUHSFklnVR8DVwGHCqyPmY25rB6VUilpJLS7/0LHtm5NuqpKpfWRqtPVzFkL4JotnOtQkYHhWhq6kSQ9V9KBytNzga9K+g7wTeBzEfGFAutjZmMu71KApv31U1Owd29yUZeS+71781VgcjIZt3j44SQK9WILz+q6haymUAEL4AqblRQRf5By7AFgZ+XxfcCLivp8M9t4qi2Clrn3snbymZzM7uY5cCD9PbXOPjvf5+S9mLfq8ipo1zevfDazsTI7m3QbnT7d0H1UK6tpsW9fdt9/nj0bjh2rT/aXuwmToVmXVzcznFpwYDCzjWd2Nrmo1nYXtbrIpr1n69b15WrHEDr5nFpZXU5Sk6jXPe/HYGbWqaxpUFLSZOlUdS+KrK6r6mh6m4ZhVpKZ2XgrIiNq7V4UaQoaV6jlwGBm1qluxxDSDGhcoVahuZLMzMZa7mlQbWg1rtAHbjGYmXUjaxpUp/smDMGGPQ4MZmZ5tHOhb7VndbNzFdE91SZ3JZmZtdK40Kx6oYf0bqNWe1Y3O1cR3VNt8nRVM7NWpqfTZwllTRttNo01azV0h1NQ2+HpqmZmvdJuzqNm4wS9yJ9UMAcGM7NW2h0QbjZOMASDy604MJiZtdLugHCzVBhDMLjcigODmVkrneZWSpvG2m3+pD7w4LOZ2QbhwWczM+uIA4OZmdXpKjBIeoOkeyWdljTT8No7JR2W9ENJr8p4/wWSviHpR5I+JWlzN/UxM7PuddtiOAT8M+ArtQclXUSy5/MLgKuBmyVNprz/z4H3RsSFwM+Bt3RZHzMz61JXgSEivh8RP0x56RrgkxHx64i4HzgMXFZbQJKA3wH+unJoH/B73dTHzMy6V9QYw3nAT2qeH60cq7UV+EVEnGxSxszM+qxlEj1JdwLPTnlpPiI+k/W2lGON82LzlKmtxxwwB7BjiFYImpmNm5aBISJe0cF5jwLn1zx/HvBAQ5mHgbMlbaq0GtLK1NZjEViEZB1DB3UyM7MciupKuh24VtJTJF0AXAh8s7ZAJCvr/g54feXQbiCrBWJmZn3S7XTV10k6CrwU+JykLwJExL3Ap4HvAV8A3hoRpyrvOSDpuZVT/DHwR5IOk4w5fKSb+piZWfdGMiWGpFUgJaH5OttIuqzGxbh9Hxi/7+TvM/zG7Tu1831KEbG9VaGRDAx5SVrKkxdkVIzb94Hx+07+PsNv3L5TEd/HKTHMzKyOA4OZmdUZ98CwOOgK9Ni4fR8Yv+/k7zP8xu079fz7jPUYg5mZtW/cWwxmZtamsQwMWenAJU1L+ntJd1duHxxkPfPqNr35sJP0bkk/rfm77Bx0nToh6erK3+GwpHcMuj7dkrQs6Z7K32Qkt0yUdKukhyQdqjl2jqQ7Kun+75D0zEHWsR0Z36fn/3/GMjCQkQ684scRcUnltqfP9epUt+nNR8F7a/4uBwZdmXZV/t0/ALwauAjYVfn7jLqXV/4mozq986Mk/zdqvQO4q5Lu/67K81HxUdZ/H+jx/5+xDAxN0oGPpG7Sm1vfXAYcjoj7IuJx4JMkfx8boIj4CvBIw+FrSNL8w4il+8/4Pj03loGhhQsk/R9J/0vSbw+6Ml3Kk958VLxN0ncrTeWRadrXGKe/RVUAX5J0sJLdeFycGxEPAlTunzXg+vRCT///jGxgkHSnpEMpt2a/0h4EdkTEbwF/BHxc0tP7U+PmOvw+baUuH6QW3+8W4DeAS0j+Rn8x0Mp2ZmT+Fm24PCJeTNI99lZJLxt0hSxVz///tEy7Paw6SQceEb8Gfl15fFDSj4F/CAx8YK3A9OZDIe/3k/Qh4H8UXJ0ijMzfIq+IeKBy/5Ck20i6y9LG7UbNzyQ9JyIelPQc4KFBV6gbEfGz6uNe/f8Z2RZDJyRtrw7OSvoHJOnA7xtsrbrSMr35KKj856x6Hclg+6j5FnChpAskbSaZFHD7gOvUMUlbJJ1VfQxcxWj+XdLcTpLmH8Yg3X8R/39GtsXQjKTXAe8HtpOkA787Il4FvAx4j6STwClgT0QUPpDTrazvExH3SqqmNz9JTXrzEfMfJF1C0vWyDPzhYKvTvog4KeltwBeBSeDWSvr5UXUucJskSK4TH4+ILwy2Su2T9AngCmCbki0CbgT+DPi0pLcAR4A3DK6G7cn4Plf0+v+PVz6bmVmdDdWVZGZmrTkwmJlZHQcGMzOr48BgZmZ1HBjMzKyOA4OZmdVxYDAzszoODGZmVuf/A8QkYe7Nga5pAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "w0=generateGalaxy()\n",
    "t= np.linspace(0,128,128)\n",
    "w = odeint(galaxyOrbit1, w0, t)\n",
    "plotGalaxy(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
